function MV = velArray(M,nhstep,nhm,vcol,str_cell)
%velArray - Computation of the velocity array
%
%       MV = velArray(M,NHSTEP,NHM,VCOL,STRCELL)
%
% Computes the velocity array MV from the position array M, with a time 
% step of NHSTEP hours and mean on NHM steps by computing the least square
% linear fit in the time span [tc-NHSTEP*NHM,tc] for each computation time
% tc, where
% 
%       MV(:,:,k) = [tc x y z vx vy vz]
%
% The 3-components vector VCOL indicates the columns of the array M to be 
% used to compute the velocity. For the k-th target, the velocity is 
% therefore computed on the basis of M(:,VCOL,k) data. 
% If VCOL is undefined or empty, the last three columns of the array M 
% (i.e. for each matrix M(:,:,k)) are considered.
% The target names are taken from the input cell variable STRCELL, which 
% must be compatible with MV, i.e. the condition size(MV,3)=numel(STRCELL)
% must be satisfied.
% If STRCELL is undefined or empty, this default cell variable is used:
%   STRCELL = {'P1','P2','P3','P4','P5','P6','P7','P8','P9','P10','P11','P12',...
%     'P13','P14','P15','P16','P17','P18','P19','P20','P21','P22','P23',...
%     'P24','P25','P26','P27','P28','P29','P30','P31',...
%     'N1','N2','N3','N4','N5','N6','N7','N8','N9','N10',...
%     'GPS_ROV 1','GPS_ROV 2',...
%     'R1','R1_BIS','R2','R2BIS','R3','R4','R5'}.
%
% See also extraTarget2cell, cell2posArray, selectMeanVelocity,
%   extraRTS2array.

warning ('off','all');

if nargin < 5
    str_cell = [];
end
if nargin < 4
    vcol = [];
end

[nt,nc,ns] = size(M);

if isempty(vcol)
    vcol = nc-2:nc;
elseif numel(vcol) ~= 3
    error('vcol mus be a 3-element vector');
else
    vcol = (vcol(:))';
end


thr = 20000; % VELOCITY THRESHOLD 

if isempty(str_cell)
    str_cell = {'P1','P2','P3','P4','P5','P6','P7','P8','P9','P10','P11','P12',...
        'P13','P14','P15','P16','P17','P18','P19','P20','P21','P22','P23',...
        'P24','P25','P26','P27','P28','P29','P30','P31',...
        'N1','N2','N3','N4','N5','N6','N7','N8','N9','N10',...
        'GPS_ROV 1','GPS_ROV 2',...
        'R1','R1_BIS','R2','R2BIS','R3','R4','R5'};
end
ns1 = numel(str_cell);   % number of stations
if ns1 ~= ns
    error('incompatible data: the number of stations is invalid');
end

dt = 1/24;          % time step of time series, in days 
dtn = dt*nhstep;    % time step of computation, in days    

if vcol(1) >= 5
    MV = nan(nt,nc,ns);
else
    MV = nan(nt,7,ns);
end
nhe = nhm*nhstep;

for k = 1:ns
    strk = str_cell{k};
    fprintf('Processing of station %s in progress\n',strk);
    Mk  = M(:,:,k);
    Ik5 = abs(Mk(:,vcol(1))) < 10^3*eps;
    Mk(Ik5,vcol(1)) = NaN;
    Ik6 = abs(Mk(:,vcol(2))) < 10^3*eps;
    Mk(Ik6,vcol(2)) = NaN;
    Ik7 = abs(Mk(:,vcol(3))) < 10^3*eps;
    Mk(Ik7,vcol(3)) = NaN;
    MVk = nan(nt,7);
    for h = 1:nt
        if h > nhe
            phm = Mk(h-nhe:h,[1 vcol]);
        else
            phm = Mk(1:h,[1 vcol]);
        end
        if size(phm,1) > 1
            ph = nanmean(phm(:,2:end),1);           % mean positions 
        else
            ph = phm(2:end);
        end
        Ivh  = ~isnan(phm(:,2));                    % valid measures
        phmr = phm(Ivh,:);
        nvh  = size(phmr,1);
        if nvh >= 3
            vh = zeros(1,3);
            for n = 1:3
                pmn = polyfit(phmr(:,1),phmr(:,vcol(n)),1); % least square straigh line
                vmn = pmn(1);
                vh(n) = vmn;
            end
        elseif nvh == 2
            vh = (phmr(2,2:end)-phmr(1,2:end))/(phmr(2,1)-phmr(1,1));
        else
            vh = nan(1,3);
        end
        Ivh = abs(vh) > thr;
        vh(Ivh) = thr*sign(vh(Ivh));
        MVk(h,2:end) = [ph vh];
    end
    MVk(:,1) = Mk(:,1);
    MV(:,:,k) = MVk(:,:);

end